Spatiotemporal Heterogeneity of Local Free Volumes in Highly Supercooled Liquid 
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We discuss the spatiotemporal behavior of local density and its relation to dynamical heterogeneity 
(DH) in a highly supercooled liquid by using molecular dynamics simulations of a binary mixture 
with different particle sizes in two dimensions. To trace voids heterogeneously existing with lower 
local densities, which move along with the structural relaxation, we employ the minimum local 
density for each particle in a time window whose width is set along with the structural relaxation 
time. Particles subject to free volumes correspond well to the mobile region of DH. However, no 
growth in the correlation length of heterogeneity in the minimum local density distribution takes 
place. 

PACS numbers: 64.70.kj, 81.05.Kf, 61.54.Fs 
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The reason for the drastic slow down in dynamics when 
liquids are cooled toward the glass transition tempera- 
ture has been a long-standing problem. In the last two 
decades, numerous efforts have been made to study this 
problem via molecular dynamics (MD) simulations [l]]. 
The recently proposed concept of "dynamic heterogene- 
ity" (DH) [2}-|6j has shown the existence of contrast be- 
tween mobile and immobile regions of supercooled liquids 
on a large scale, as if there were critical fluctuations hid- 
den behind. In recent literature, considerable attention 
has been focused on relation of DH to the structural het- 
erogeneity of medium length-scale crystalline order 
icosahedral order local potential energy (lfjj. and so 
on. In the context of experiments especially on metallic, 
polymer, and colloidal glasses, density fiucstuations have 
been observed over short to long lengths [llHl3|. One 
of the primary theoretical challenges had been the con- 
struction of statistical descriptions of solids and glasses 
including the effects of vacancies or local- free volume[l4r- 
Il7j . In MD simulations of binary mixtures usually em- 
ployed in studies of supercooled liquids, density fluctua- 
tions are too weak to be captured clearly 0, U 18 1, and 
the correlation of such fluctuations to DH has only par- 
tially been observed flij]. In two-dimensional (2D) crys- 
tals under melting, where the density fluctuation is more 
pronounced than in glass, the diffusion of defect clusters 
with relatively smaller local densities is ascribed as the 
cause of DH|20| . In this Letter, we clarify the relationship 
between the local density and DH by capturing and trac- 
ing weakly existing density fluctuations in a constructive 
manner. 

We investigate a 2D binary mixture composed of two 
atomic species, 1 and 2, with N\ = N 2 — 32000 parti- 
cles. The particles interact via the soft-core potentials 
v a p{r) = e(cr Q(3 /r) 12 + C Q/3 where cr Q(3 = (a a + crp)/2 
and r denote the interaction lengths and the distance 
between two particles respectively, with a, j3 € {1,2}. 



The interaction is truncated at r = 4.5cri and a constant 
C a p is set so as to ensure continuity of the potential 
at the cutoff. The size ratio between the two species 
is <J<il<T\ = 1.4 to prevent crystallization, and the mass 
is set as m<xjm\ = (02/eri) 2 . The particle density is 
fixed at a high value of <j> = (Niaf + N 2 cr%)/V = 1.2, 
and therefore the particle configurations are jammed in 
the supercooled state. No tendency of phase separation 
is detected in our computation times. Space, t ime, and 
temperature are measured in units of a%, To = y niiaf/e, 
and e/ks, respectively. Dependence of a- relaxation time 
on the temperature T is shown in Table HI A sufficiently 
long annealing time tA is chosen (tA > 30r a ) with the 
time step being At = 0.005. No aging effect was ap- 
preciable in the course of calculation of pressure, density 
time correlations, etc. After performing the above pro- 
cedure at each temperature, we begin to collect the data. 
This time point in each of our simulations is denoted as 
t = in the following. 

First, for each particle i, we define the local density by 
counting the particles within the distance Ar = \r — r%\ < 
ro from i 



M = zi2 



dr <Jj5(rj 

Ar<ro 
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This local density is weighted by the assumed area cr|. 
The inset of Fig. 1 shows the long-time average of the 
local density [n,(ro)]iv at T = 0.56. Here, [AJw = 
TV -1 J^i Ai denotes the average over all the particles. 
This quantity corresponds to a radial distribution func- 
tion g(ro) defined in a naive manner. 

For the same simulation run, at appropriate time inter- 
vals, we estimate the local density in the mobile regions of 
DH as follows: bonds are defined at each time to as par- 
ticle pairs between i and j satisfying rij(io) < Aia a p, 
and after a time interval At bonds are regarded to be 
broken if r^-(to + At) > A2<J a p [2], 0]. The cutoffs are 
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TABLE I: Dependence of a-relaxation time (r a ) and bond-relation time (17,) on T. 



T 


0.56 


0.64 


0.72 


0.80 


0.96 


1.20 




2.14 x 10 3 


75.5 


10.5 


5.39 


2.86 


1.47 


n 


1.41 x 10 5 


1.90 x 10 4 


5.25 x 10 3 


1.87 x 10 3 


5.35 x 10 2 


1.85 x 10 2 




FIG. 1: Degree of deviation R(ro) of the local densities around 
the broken bonds [w(ro)]b.b. (obtained at intervals of At — 
2 x 10 2 ) from that averaged over the entire system [rij(ro)]iv, 
plotted as a function of averaging radius ro. The data for 
T — 0.56,0.64, and 0.72 are shown. Inset: (Ui(ro)) is shown 
for T = 0.56. 




FIG. 2: The red points show the distribution of 2% of the 
particles having smaller local densities {rii(ro)}m, which is 
averaged over t € [0, 10], in an IE of 32 runs at T = 0.56. 
The number distribution of broken bonds for the same IE 
(B(r, At)) ie with t € [0, 122] is shown according to the scale 
bar on the right. For clarity, a magnified figure is shown on 
the bottom right. 



set to A\ = 1.15 and A 2 = 1.6, respectively. We take 
the local density profile only for particles that have un- 
dergone these broken bonds for At = 2 x 10 2 , and we 
denote their average local density profile as [nj(»"o)]b.b.- 
The main graph in Fig. 1 shows the long-time average 
R(ro, At) of the relative degree of deviation as given by 
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FIG. 3: The red solid line shows the local density rii(ro) (ro = 
6.0ai ) averaged over every unit time (Atd = 1) for a run 
at T = 0.56. On the blue dashed line, ni(ro) is calculated 
for every time interval of At^ = 0.05 and each minimum is 
plotted per unit time. Therefore, the latter is systematically 
lower than the former due to thermal fluctuations on a short 
time scale. Inset: For the same particle and for a longer time 
interval, m(ro) is plotted when averaged over a time interval 
of At d = 10. 



While R(ro) exhibits an oscillating behavior at small val- 
ues of ro mainly due to the mixed contributions from 
larger and smaller components, its value is systemati- 
cally lower than zero when ro is sufficiently large. Thus, 
around the broken bonds, i.e., in the mobile regions, the 
smaller local density can be characterized with a suffi- 
ciently large value of ro. We employ ro = 6.0 in our 
following discussions. 

To characterize the mobile regions of DH, we count 
the number of broken bonds pairs for each particle as 
Bi(At) = X)jeb b 1; where the summation is taken over 
j particles at the broken bond ends of i particles cal- 
culated over time intervals of length At. Since DH is 



R(r ,At) = 



N(r )]b.b. 
[ni(r )] N 



- 1 



(2) 



uniquely dependent on the particle positions 2l| , isocon- 
figurational ensembles (IEs) of 32 runs are employed and 
an average is obtained over these runs. This IE average 
is denoted as (-)ie in the following discussions. In Fig. 
[H we demonstrate (S,(At)) IB for a short time interval, 
t G [0, At] (At = 122), together with distribution of 2% 
of particles having the lowest local densities (rij(ro))iE 
taken at the initial stage t € [0, 10]. Clusters of particles 
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FIG. 4: The red points show the distribution of 25% of the particles having smaller values of minimum local densities {vi(At b ))iE, 
in an IE of 32 runs performed for Atb = 6,60, and 6100 for T — 1.2,0.80, and 0.56, respectively. The number distribution of 
broken bonds (Bi(Atb))m is shown by the gray scale according to the scale bar on the right. 



with low local densities are observed, and the particles 
are clustered at distances of several times the particle ra- 
dius. On a mesoscopic scale, this clustering corresponds 
to particle propensity represented by (B(r, At))iE, with 
B(r,At) = J2j^jBj(At)5{r - rj(0)). From this corre- 
spondence, a natural speculation would be that the local 
free volumes within each cluster intract with each other 
with a blurred heterogeneity and move around, thereby 
leading to a clear heterogeneity in the dynamics. 

To understand the cluster relationship with DH, it is 
necessary to capture how they move over longer time 
intervals. For T = 0.56,0.64,0.72,0.80,0.96, and 1.20, 
simulation runs of IEs arc performed with time inter- 
vals of At b (T) = 6100, 815, 210, 80, 23, and 8, re- 
spectively. Upon using the concept of bond relaxation 
time Tfe(T) introduced in a previous study [H and defined 
by Nb(to + n) = N b (t )/e, these time intervals satisfy 
At b (T) ~ 0.043t 6 (T). Here, N b (t) denotes the total num- 
ber of bonds at time t, and T b {T) characterizes the time 
scale of structural relaxation (see Table U for actual val- 
ues). In these intervals, almost identical portions of the 
total initial bonds become broken for all values of T. 

At low T where slow dynamics realize, particle rear- 
rangements agitated by the thermal fluctuations are ex- 
pected to become intermittent. In Fig. [H for T = 0.56, 
change in local density rz.j with time is measured for one 
specific particle in the mobile region of DH. Averaging n% 
over every time interval of Atd = 10 as shown in the in- 
set, the local density is observed to be fluctuating around 
4>. In the region around t = 2700, the density reduces 
for a time interval of order 10 2 by an amount of around 
An, ~ 0.05. Interestingly, since Am x 7r(r /2) 2 ~ O(l), 
this reduction corresponds to a free volume of one parti- 
cle size scale for the current radius rg = 6.0. The lower 
the temperature, the longer the time interval is expected 
to be due to smaller thermal fluctuations. Upon decreas- 
ing the averaging time Atd to 1 and 0.05, the temporal 
fluctuation of nj(ro) appears more explicitly In the main 



graph, the red solid line represents nj(ro) averaged over 
Atd = 1- For a much shorter time interval of Atd = 0.05, 
we accumulate the data of the density, and for every in- 
terval of unit time length (At = 1), we take the minimum 
of n, (ro). The blue dashed line represents the minimum 
value of the local density for each unit time thus taken. 
Because these two lines behave in a parallel manner when 
the density drops down, the minimum value in the local 
density history assumed by one particle can well repre- 
sent and characterize the general time-development of the 
local density. 

To capture the trace of free volumes and its correspon- 
dence to DH, we investigate the spatial distribution of 
the "minimum local density" that is defined as 

vAAt) = min rii, t\ = to + At. (3) 

te[t ,ti] 

The above term indicates the minimum value in the 
density history for each particle i. The IE average of 
u(r,At b ) = J2jVjS{r — Tj(0)) can illustrate the spa- 
tiotemporal heterogeneity of free volumes in an enhanced 
manner. Figure 0] shows the distribution of 25% of the to- 
tal number of particles having lower (i/j (At b )}iE mapped 
together with the broken bond distribution (Bi(At b ))m 
The lower the value of T is, the more distinct is the cor- 
respondence between these distributions. 

Quantitatively, this correspondence is illustrated for 
T = 0.56 in Fig. [5] the correspondence is observed as 
a scattered plot between (0(r, At b ))\B, with u(r,At b ) = 
J2j VjSlr — rjfi)) and (B(r, At b ))iE averaged in each cell, 
thereby squarely dividing the total system into a 12 x 12 
grid of size (L/12) 2 . 

To examine how the heterogeneity in the density 
evolves, in Fig. [51 we show the structure factor of "free 
volume degree" 0(r, At b ) defined as 

S v (k,At b ) = (\0 k (At b )\ 2 ) m , (4) 

where 80(r,At b ) = V ; rr;,//, - [iSj} N )5(r - r\,-(0)) is 
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FIG. 5: Scattered plot of cell-averaged quantities between 
{P(r, At b )) m and (B(r, At b ))i E for T = 0.56. The cells divide 
the total system into a 12 x 12 grid. The data for runs with 
eight independent initial conditions are shown in this figure, 
and an IE average over 32 runs is taken for each run. 




FIG. 6: Structure factor of the "free volume degree" 
S v (k, Atb) defined in Eq. @ for various values of T. Inset: 
Structure factor of the broken bonds S b (k, At b ), exhibiting 
the Ornstein-Zernike form typical to DH whose correlation 
length grows as T decreases. 

the mesoscopic fluctuation of Vi, and i>k{Atb) indicates 
Fourier components of Si)(r, At b ). The average is taken 
over IEs of 32 runs generated from 32 independent initial 
configurations, and thus, data from 1024 runs performed 
over a duration of Atb are used for each T. In the inset of 
Fig. |B1 we also show the structure factor of corresponding 
broken bonds 

S b (k,At b ) = (\B k {At b )\ 2 )iE (5) 

which gives information on how DH evolves at long- 
wavelengths. Since one isolated free volume lower the val- 
ues of rii(ro) of particles wthin the distance of I = tqGi, 
we have a small peak of S„(k,Atb) around k ~ 27r// ~ 
10°, and this peak is enhanced at long- wavelengths. Re- 
markably, S v (k, At b ) does not change with tempera- 



ture exhibiting the same degree of heterogeneity even at 
long- wavelengths, while Sb[k, At b ) displays the Ornstein- 
Zernike form with growing length scales for lower values 
of T in a manner similar to those obtained in previous 
studies[3, Q. That is, for the interval of corresponding 
structural relaxation time At b (T) the heterogeneity in 
i>(r, At b ) is nearly independent of the temperature, and 
B(r, Atb) shown in Fig. |4] might suggest existence of 
other structural origins for DH; as regards the arguments 
on these origins in the literature, such as density gradi- 
enty, growing static length scales, and so on, their rela- 
tion to our results should be examined in the future. 

In conclusion, we have clarified the relation between 
free volume distribution and DH for various temperature 
for the first time, using MD simulations of a binary mix- 
ture. While the distribution of free volumes is extended 
at a mesoscale to form clusters, within each cluster, sev- 
eral free volumes are scattered at the one-particle level, 
as shown in Fig. [2l These complex characteristics have 
thus far not been discussed in the literature of MD sim- 
ulations of glasses. DH can be attributed to cooperative 
rearrangements of several free volumes in such clusters 
due to the effect of thermal excitations. In our binary 
mixture, the density distribution and the relative degrees 
of local density drop down at a free volume are suggested 
to be independent of the temperature, with no apparent 
length scales growing with decreasing temperature. 
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